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Chapter  1 

Surrogate  models 


1.1  Introduction 

Design  of  complex  microwave  systems  is  carried  out  using  commercial  circuit  simulators.  This 
is  a  long  process  because  the  design  involves  the  optimization.  In  microwaves  circuits  are 
distributed,  so  the  response  depends  on  the  dimensions  and  topology.  The  response  can  be  eval¬ 
uated  using  either  the  closed  form  formulas  or  full-wave  electromagnetic  simulators.  Full-wave 
simulators  are  too  slow  to  be  used  in  an  optimization  loop,  so  all  circuit  simulators  employ 
closed  form  formulas,  or  in  other  words  surrogate  models.  The  accuracy  and  speed  of  anal¬ 
ysis  of  a  simulator  relies  on  fast  and  accurate  surrogate  models  of  microwave  discontinuities. 
Surrogates  are  used  on  many  stages  of  a  microwave  circuit  design:  from  initial  design  to  final 
optimization  and  yield  analysis.  For  a  microwave  designer  one  of  the  most  import  requirements 
regarding  the  circuit  simulator  is  the  quality  and  diversity  of  models  library.  Recently,  new  tech¬ 
nologies  have  emerged  for  low  cost  millimeter-wave  systems  based  on  the  system  on  the  chip 
(SoC)  or  system  on  the  package  (SoP)  design  philosophies.  For  instance,  SoC  employs  silicon 
germanium  (SiGe)  on  either  CMOS/BiCMOS-grade  Si  or  highresistivity  Si  as  a  replacement 
for  GaAs  for  some  applications.  SoP  modules  integrate  different  passive  components  in  a  mul¬ 
tilayer  low  loss  material  such  as  low  temperature  co-fired  ceramic  (LTCC)  or  Liquid  Crystal 
Polymers.  Other  category  are  millimeter  wave  devices  made  in  multilayer  MCM-D  technology. 
For  these  emergent  technologies  there  are  no  component  libraries  that  can  be  used  with  the  in¬ 
dustry  standard  circuit  simulators  such  as  Agilient’s  ADS  or  AWR’s  Microwave  Office.  As  a 
result,  one  has  either  to  apply  approximate  formulas  that  are  often  too  inaccurate  or  resort  to 
a  full  wave  solver  in  order  to  design  systems  components.  EM  simulations  are  too  computa¬ 
tionally  intensive  to  be  used  in  the  optimization  loop.  The  design  time  would  be  significantly 
reduced  if  a  designer  could  use  the  surrogate  parameterized  models  that  could  be  evaluated 
at  the  speed  of  closed  form  formulas  but  having  the  accuracy  comparable  to  EM  simulations. 
One  of  the  goals  of  the  research  in  this  area  is  to  create  a  technique  that  allow  one  to  build  a 
multidimensional  parametric  model  of  a  component  using  as  few  data  points  as  possible. 

The  development  of  such  a  technique  was  the  main  motivation  for  this  research.  The  ba¬ 
sic  algorithm  of  surrogate  model  construction  was  developed  previously  and  published  in  two 
publications  [1,  2].  The  main  assumption  is  to  represent  the  transfer  function  of  the  device  be¬ 
ing  modelled  with  a  multivariate  rational  function  with  adaptive  support  point  and  model  order 
selection.  The  technique  is  the  extension  of  the  technique  presented  in  [15]  to  the  multivariate 
case.  The  adaptive  sampling  over  whole  parameter  space  was  utilized  to  efficiently  select  and 
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limit  the  number  of  support  points.  Detailed  description  of  the  technique  will  be  given  in  next 
sections. 

Given  the  bacground  and  the  fact  that  the  basic  algorithm  was  already  available,  the  goals 
of  this  project  was  to  advance  the  published  technique  and  prove  its  viability  for  microwave 
circuit  design  by  creating  multivariate  surrogate  models  of  high  complexity  of  componetns  in 
one  of  the  emerging  technologies  and  ready  to  use  in  industry  standard  circuit  simulators  (ADS, 
Microwave  Office).  The  models  have  the  accuracy  comparable  to  full  wave  simulations  but  at 
the  same  time  the  computational  speed  similar  to  the  closed  form  formulae.  As  a  result  one  is 
able  to  achieve  fast  optimization  of  microwave  circuits  manufactured  in  the  emerging  and  new 
technologies.  Specifically  the  following  issues  were  addressed: 

•  Improvement  of  stability  of  the  interpolation  solvers  by  replacing  multinomials  with  bet¬ 
ter  conditioned  orthogonal  polynomials. 

•  Development  of  the  technique  for  automated  division  of  the  parameter  space  that  will 
allow  on  to  create  different  low  order  models  in  each  subspace 

•  Development  of  the  technique  for  merging  the  submodels  into  a  single  model  covering  a 
wider  parameters  range. 

•  Integration  of  the  software  for  automated  model  generation  with  the  industry  standard 
planar  structure  simulators  such  as  e.g.  Sonnet. 

•  Development  of  software  for  automated  generation  of  compiled  models  that  can  the  used 
with  ADS  and  Microwave  Office 

•  Demonstration  of  the  suitability  of  the  approach  for  emerging  technologies  such  as  LTCC 
and  LCP,  MCM-D  by  developing  models  of  typical  discontinuities  or  elements  e.g.  spiral 
inductors. 


1.2  Alternative  solutions  for  model  construction 

The  technique  chosen  for  this  effort  as  a  basis  for  the  surrogate  model  construction  has  a  few 
alternatives.  One  popular  solution  which  is  applied  for  surrogate  models  involves  artificial 
neural  networks  [3,  4,  5],  however  the  drawbacks  of  ANN  (unknown  network  topology  and 
long  training  process)  significantly  limit  their  usage  in  automated  model  construction.  The 
surrogates  can  also  be  created  with  application  of  radial  basis  functions  (RBF)  [6,  7].  The  issue 
is  the  selection  of  best  value  of  unknown  shape-parameter  of  radial  functions  [8].  On  the  other 
hand  the  approach  using  RBF’s  significantly  reduces  the  problems  with  ill-conditioning.  Test 
carried  out  by  our  research  group  showed  however  that  the  RBF’s  are  inefficient  in  case  of 
complex  devices. 

Another  approach  for  automatic  model  creation  was  presented  in  [9].  In  this  algorithm,  fre¬ 
quency  is  handled  separately  from  other  physical  parameters.  The  procedure  has  two  stages:  at 
first  at  selected  frequency  points  multidimensional  models  are  created  by  expanding  the  multi¬ 
variate  functions  into  series  of  orthogonal  multinomials.  The  expansion  coefficients  are  found 
by  solving  a  system  of  interpolator  conditions.  On  this  stage  the  support  points  are  added  in 
an  entirely  adaptive  way.  Next  the  frequency  dependence  is  added  by  one-dimensional  rational 
interpolation  of  the  models  response.  The  procedure  creates  models  with  good  accuracy,  but  it 
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is  obvious  that  excluding  frequency  from  the  adaptive  sampling  procedure  may  result  in  non- 
optimal  number  of  support  points.  Presented  results  also  show  that  such  solution  is  efficient  up 
to  three  parameters. 

Lehmensiek  and  Meyer  [10,  11,  12]  developed  techniques  based  on  Thiele-type  branched 
continued  fraction  representation  of  a  rational  function.  The  algorithms  operate  by  using  a  uni¬ 
variate  adaptive  sampling  along  a  selected  dimension.  In  this  way,  while  the  support  points  do 
not  fill  the  grid  completely,  they  are  being  added  along  straight  lines  passing  through  multidi¬ 
mensional  space.  The  efficiency  of  the  algorithms  was  illustrated  on  two-  and  three-dimensional 
models. 

Yet  another  approach  uses  statistical  tools  for  model  construction  like  Kriging  [13]  and  De¬ 
sign  of  Experiment  (DOE)  [14]  techniques.  Kriging  is  a  special  form  of  interpolation  function 
that  employs  the  correlation  between  neighboring  points  to  determine  the  overall  function  at  an 
arbitrary  point.  DOE  makes  a  series  of  tests  in  which  a  set  of  input  variables  is  changed  and 
the  experimenter  can  identify  the  reasons  for  changes  in  the  output  response.  Based  on  this 
knowledge  one  can  construct  a  statistical  model  if  the  test  structure.  Both  techniques  can  be 
applied  to  create  of  simple  models  with  low  accuracy  and  are  dedicated  to  coarse  tuning  of  the 
design. 


1.3  Measures  of  model  quality 

In  every  model  construction  scheme,  the  basic  question  is  how  to  assess  the  accuracy  of  the 
created  model.  Estimation  of  the  model  error  is  an  important  issue,  because  in  many  cases  the 
accuracy  of  model  limits  the  range  of  its  applications.  There  are  several  criterions  of  verification 
of  model  accuracy.  In  techniques  that  involves  adaptive  sampling,  like  one  described  in  this 
report,  the  basic  measure  of  error  is  the  maximum  difference  between  two  different  models 
which  are  used  to  select  the  set  of  support  points.  Let  us  define  this  error  as  £.  As  described 
later,  the  goal  of  the  procedure  is  to  create  two  different  models,  in  case  of  which  the  error  £  is 
below  required  accuracy  £o.  However  the  criterion  £  <  £o  does  not  guarantee  that  created  model 
have  accuracy  as  good  as  £o-  To  estimate  the  real  accuracy  of  the  model,  one  needs  to  perform 
an  additional  statistical  test.  The  test  involves  computation  of  model  response  on  a  test  set  of 
random  points  scattered  in  the  model  parameter  space  and  compare  the  results  with  response  of 
device  being  modelled.  In  result  one  gets  a  real  accuracy  of  the  model,  defined  in  this  report  as 
A. 

Form  the  data  evaluated  and  randomly  selected  points  one  can  compute  various  statistical 
measures.  In  this  report  we  use  the  following  simple  ones:  the  maximum  real  error  Amax  and 
mean  real  error  Amean.  Both  errors  can  be  computed  in  decibels  as  A[dB]  =  201og(A).  Such 
maximum  and  mean  real  errors  give  a  good  measure  of  model  accuracy.  They  are  however  too 
simple  for  practical  purpose,  because  they  do  not  provide  the  information  regarding  the  expected 
accuracy  for  a  particular  set  of  parameters.  Therefore,  it  appears  that  it  is  more  meaningful  to 
apply  the  cumulative  distribution  function  of  error  A  to  derive  the  quality  measure.  In  our  case 
we  descided  to  use  the  error  level  A90  that  fulfills  condition  that  90%  of  points  in  the  test  set  have 
an  accuracy  better  than  A90.  As  discussed  in  [16]  the  error  estimate  £  has  a  better  correlation 
with  A90  than  with  Amean  or  Amax. 
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1.4  Technique  details 


The  technique  selected  for  this  work  creates  an  interpolant  of  N  variate,  real  or  complex  valuated 
smooth  function  S (X)  =  S(xi  ,*2,  •  •  • ,  xn)  as  a  rational  function  [  1  ] : 


S(x  i,x2,...,xN) 


AQQ  _  A{x i,x2,...,xN) 
B(X)  B(xi,x2,...,xn) 


(1.1) 


where  both  numerator  A (X)  and  denominator  B(X_)  are  multinomials  (sum  of  monomials  mul¬ 
tiplied  with  scalar  coefficients).  The  complete  set  of  the  monomials  can  be  listed  as  elements  of 
matrix[17]: 


Row  1  :  1 

Row  2  :  x\  X2  ...  xn 

Row  3  :  xj  X1X2  . . .  x  1  x,\>  x\  X2X3  . . .  x\. 

Row  4  :  x\  x\X2  v^.v'3  . . .  x\x^  x\  x\x  \  x\x3  . . .  xjj 


In  the  m-th  row  are  written  all  monomials  with  sum  of  powers  at  each  variable  equal  m-1 .  With 
such  assumption  the  multinomials  of  numerator/denominator  can  be  described  with  a  vector 
V  =  [v [ ,  v’2,  •  ■  ■ ;  v’/v],  where  v,-  determines  maximum  allowed  power  of  the  i-th  variable.  For 
example  in  case  of  a  three- variate  multinomial  which  order  is  described  with  vector  V  —  [3  2  2] 
the  following  monomials  are  selected: 

Row  1  :  1 

Row  2  :  x\  X2  X3 

Row  3  :  x\  x  1 V2  X1X3  x\  X2X3  x^ 

Row  4  :  x\  X\X2  x\x3  X\x\  x\x3  X 1 V3  X2X3  X 1X2X3 

In  further  investigations  it  is  assumed  that  both  numerator  A  and  denominator  B  of  (1.1)  have 
the  same  orders,  therefore  Va  —  Vb  —  V. 

The  unknown  coefficients  a,-  and  b,  corresponding  to  the  multinomials  of  numerator  and 
denominator  of  1.1  can  be  found  requiring  that  equation: 

A(X)-5(X)R(X)  =0  (1.2) 


is  fulfilled  in  at  least  L>M\  +M2  support  points,  where  M 1  and  M2  are  the  numbers  of  unknown 
coefficients  a,-  and  bt.  The  fitting  problem  can  be  transformed  to  the  matrix  form: 


[A  -B] 


(13) 


where  a  and  b  are  the  vectors  of  unknown  coefficients  and  [A\lxMx  ,  [B]lxM2  are  matrices  involv¬ 
ing  the  values  of  the  monomials  appearing  in  numerator  and  denominator  of  (1 . 1)  as  well  as  the 
values  of  the  interpolated  function  at  support  points. 

The  problem  can  be  solved  applying  the  total  least  squares  technique,  as  described  in  [18]. 
The  total  least  squares  method  is  suited  to  problems  in  which  both  the  coefficient  matrix  and 
the  right-hand  side  are  not  precisely  known.  It  allows  one  to  filter  the  noise  from  interpolated 
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data  and  improve  the  quality  of  resulting  solution.  The  first  step  is  the  computation  of  the  QR 
decomposition  of  the  matrix  C  =  [A  —B\,  which  results  in: 


^11  ^12 

a 

0  R22 

b 

-0 


(1.4) 


where  matrix  7?n  has  size  (M\  x  M\),  R\2  has  size  (. M\  x  M2)  and  R22  is  (( L  —  M\  +  1)  x  M2) 
matrix.  The  R22  matrix  is  affected  by  the  noise.  The  equation  (1.4)  can  be  written  as  two 
separate  equations: 

[Rn]a  =  -[Rl2\b  (1.5) 

[R2l\b  =  0  (1.6) 


Computing  the  singular  value  decomposition  ( SVD )  of  R 22  one  obtains: 


[[/][£]  [1^  =  0 


(1.7) 


where  matrix  V  is  size  (M 2  x  M2).  The  solution  of  the  problem  in  TLS  sense  is  proportional  to 
the  last  column  of  the  matrix  V,  therefore  is  assumed  that  b  —  [V] ,.v/2- 


1.4.1  Condition  number  improvement 

Condition  number  in  least  squares  measures  the  sensitivity  of  the  solution  of  a  system  of  linear 
equations  to  errors  in  the  data.  The  value  of  condition  number  allows  one  to  decide  if  the 
solution  of  the  least  squares  is  reliable  and  accurate.  Namely  a  value  of  the  conditioning  number 
near  1  indicates  well-conditioned  least  squares  problem.  The  conditioning  number  is  computed 
as  a  ratio  of  the  largest  singular  value  of  matrix  that  forms  a  least  squares  problem  to  the  smallest 
one. 

A  major  issue  of  rational  interpolation  is  a  poor  conditioning  of  equation  system.  In  order 
to  cope  with  this  this  problem  two  techniques  are  recommended: 

•  Mapping  each  of  variables  to  a  line  segment  <-1,1  > 


•  Substitution  of  simple  monomials  with  orthogonal  Tchebychev  polynomials 


The  first  technique  is  linear  mapping  of  model’s  domain  to  the  multidimensional  box  with  side 
of  line  segment.  The  mapping  of  a  z'-th  variable  x ,•  is  expressed  by  formula: 


Xj  /n  2 


(Xj-XQ ,;) 

Axj 


(1.8) 


where  x,-,m  is  a  mapped  variable,  xq ,/  denotes  center  point  of  the  parameter  range  and  Ar,  denotes 
the  width  of  the  parameter  range. 

To  improve  the  conditioning  of  the  interpolation  problem  the  regular  elements  x "  composing 
the  monomials  are  replaced  with  Tchebychev  polynomials  Tn(x[),  that  are  orthogonal  on  line 
segment  <  —  1, 1  >.  Tchebychev  polynomials  can  be  computed  using  recurrence  formula: 


Tl(x) 

=  1 

Tl  M 

=  X 

O 

=  2x  —  1 

Ti(x) 

=  4x2-3x 

T„+ i(x) 

=  2x  ■  Tn(x)  -  Tn_i  (x 
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Figure  1.1:  Iris  in  rectangular  waveguide 
Table  1.1:  Range  of  parameters  of  rectangular  iris  model 


Parameter 

Range 

frequency  / 

11.855GHz-  18.02GHz 

width  a 

6.32mm  -  15.8mm 

height  b 

4.74mm  -  7.899mm 

thickness  d 

0.2mm  -  2mm 

As  an  result,  the  set  of  monomials  is  transformed  into  form: 

Row  1  :  1 

Row  2:  7i(*i)  Ti(x2)  ...  7)  (xN) 

Row  3  :  T2{x\)  T\(x\)T\(x2)  . . .  Tx{x\)Tx{xN)  T2(x2)  Ti(x2)Ti(x3)  . . .  T2 (xN) 
Row 4:  T3(x i)  T2(xi)Ti(x2)  T2(xi)Ti(x3)  ...  T2{x\)T\{xN)  T3(x2)  T2(x2)T\(x\) 
T2{x2)T\{x3)  ...T3(xn) 


Example.  The  improvement  of  the  conditioning  of  the  problem  in  case  of  a  rational  model 
of  scattering  parameter  Sn  of  an  iris  (figure  1.1)  in  rectangular  waveguide  WR62  is  presented. 
The  selection  of  the  structure  is  motivated  by  fast  computation  of  electromagnetic  response  with 
a  mode-matching  technique  and  complex  (resonant)  response  (figure  1.2).  The  model  involves 
four  parameters:  frequency,  iris  width,  height  and  thickness.  The  range  of  parameters  is  pre¬ 
sented  in  table  1.1.  The  data  were  interpolated  with  orthogonal  and  non-orthogonal  functions, 
for  different  rectangular  grid  resolution  D  and  with  or  without  mapping  of  model  domain.  The 
results  of  such  conditioning  test  are  presented  in  tables  1.2  and  1.3.  It  can  be  seen  that  the  pro¬ 
posed  approach  provides  a  significant  reduction  of  the  condition  number,  which  assures  better 
reliability  of  constructed  models. 

1.4.2  Selection  of  support  points 

Optimal  support  point  selection  is  an  essential  issue  of  every  interpolation  scheme.  It  is  espe¬ 
cially  important  in  case  of  modelling  of  multivariate  functions  where  the  number  of  support 
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Figure  1.2:  Sample  Sn(f,a)  response  of  iris  in  WR62  rectangular  waveguide  with  iris  height 
b  =  6.32mm  and  iris  thickness  d  —  1.1mm. 


Table  1.2:  Condition  number  computed  in  waveguide  iris  case  for  rectangular  grid  with  divi¬ 
sions  D  —  4  and  different  polynomials 


Model  order  V 

Space,  Polynomials 

Non-mapped,  Regular 

Mapped,  Regular 

Mapped,  Tchebychev 

[22  22] 

6692.3 

186.3 

136.5 

[3  3  3  3] 

508300 

1489.5 

1139.1 

Table  1.3:  Condition  number  computed  in  waveguide  iris  case  for  rectangular  grid  with  divi¬ 
sions  D  =  6  and  different  polynomials 


Model  order  V 

Space,  Polynomials 

Non-mapped,  Regular 

Mapped,  Regular 

Mapped,  Tchebychev 

[2  2  2  2] 

7428 

194.1 

142.7 

[3  3  3  3] 

275350 

1037 

746.4 

points  can  be  enormous.  Since  each  point  corresponds  to  one  electromagnetic  simulation  of 
device  being  modelled,  it  is  obvious  that  a  minimization  of  samples  number  is  critical.  The 
simplest  choice  of  dense  multidimensional  rectangular  grid  is  not  optimal,  due  to  fast  growth  of 
number  of  samples  with  increase  of  grid  resolution  and  number  of  model  parameters  (see  table 
1.4). 

In  proposed  technique  an  adaptive  sampling  technique  called  also  reflective  extrapolation  is 
used  [9].  The  idea  of  adaptive  sampling  is  to  create  two  different  models  (.S'  i  (X_)  and  53  QQ) 
of  the  modelled  function  and  place  additional  samples  at  the  points  of  biggest  mismatch  be¬ 
tween  both  models  £  =  ||52(X)  —  5 1  (X)  || .  Such  a  procedure,  reiterated,  leads  to  model  quality 
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Figure  1.3:  Two  dimensional  example  of  support  points  selection  -  added  points  are  placed  in 
non-regular  manner  over  the  parameter  space. 


improvement  and  minimizes  the  number  of  samples  used,  which  is  especially  advisable  if  the 
models  are  based  on  results  of  electromagnetic  simulations.  The  procedure  is  reiterated  until 
the  mismatch  £  between  low  and  high  order  model  drops  below  required  tolerance  £o- 

An  important  issue  is  a  search  procedure  of  the  point  of  biggest  mismatch.  This  problem  cor¬ 
responds  to  problem  of  search  of  maximum  of  multivariate  function  over  a  multi-dimensional 
box  and  especially  in  the  case  of  high  number  of  model  parameters  this  issue  has  to  be  solved 
efficiently.  To  deal  with  this  problem  a  genetic  optimization  procedure  was  applied.  The  ad¬ 
vantage  of  genetic  optimization  is  that  it  allows  to  find  a  global,  not  local,  minimum  of  the 
function.  The  search  is  performed  over  the  whole  parameter  space  of  the  model,  like  presented 
in  figure  1.3,  and  the  technique  gives  better  results  than  approaches  in  which  the  parameter 
space  is  searched  with  structured  pattern  (like  grids). 

It  has  to  be  noted  that  since  a  goal  of  adaptive  sampling  is  the  minimization  of  error 
£  =  1 1 5*2 (X)  —  5|  QQ||,  the  overall  accuracy  of  resulting  models  can  be  worse  than  £q.  Such 


Table  1.4:  Number  of  support  points  needed  for  V-dimensional  rectangular  grid  with  resolution 
D 


Variables 

Divisions 

1 

2 

3 

4 

5 

6 

1 

1 

2 

3 

4 

5 

6 

2 

2 

4 

9 

16 

25 

36 

3 

3 

8 

27 

64 

125 

216 

4 

4 

16 

81 

256 

725 

2196 

5 

5 

32 

243 

1024 

3125 

7776 

6 

6 

64 

729 

4049 

15625 

46656 

7 

7 

128 

2187 

16384 

78125 

279936 
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Figure  1 .4:  Range  of  model  error  for  mismatch  between  models  equal  error  £. 


situation  is  possible  when  both  models  converge  to  the  function  which  slightly  differs  from 
electromagnetic  response.  Let  A$i  and  A52  represent  the  real  error  of  models  5i(X)  and  S2GO 
related  to  electromagnetic  response.  Assume  that  the  maximum  acceptable  error  between  both 
models  is  £o-  In  this  case  the  possible  values  of  model  error  A51  and  A52  are  illustrated  in  figure 

1.4  (dotted  area).  For  example,  in  case  of  the  point  A,  the  real  error  A51,  A52  of  both  models  is 
higher  than  £,  however  the  relative  error  between  both  models  is  below  £  because  the  real  errors 
of  both  models  have  the  same  signs. 

During  adaptive  sampling  process  a  clustering  of  support  points  can  occur,  i.e.  the  subse¬ 
quent  support  points  are  added  in  the  same  location.  The  interpolation  problem  is  then  expanded 
with  the  same  points  which  do  not  give  any  extra  information  about  the  device  response,  but 
make  the  problem  bigger  and  harder  to  solve.  Such  situation  should  be  eliminated,  therefore  if 
algorithm  detects  such  behavior,  the  parameter  space  is  divided  into  2N  smaller  subspaces  (each 
dimension  is  divided  to  half)  and  the  locations  of  biggest  mismatch  between  models  in  those 
2n  subspaces  is  found.  The  set  of  points  is  appended  to  the  data  set  and  the  adaptive  sampling 
continues. 

Implementation  of  described  adaptive  sampling  makes  it  possible  to  detect  if  the  interpola¬ 
tion  problem  is  ill-conditioned.  The  models  obtained  as  a  solution  of  ill-conditioned  system  do 
not  match  each  other  and  in  result  error  between  both  models  is  large  (namely  £  >  1).  If  such 
situation  occurs  on  the  initial  stage  of  model  construction,  when  the  number  of  support  points 
is  small,  the  solution  is  to  add  more  points  to  the  system  and  improve  the  conditioning.  Another 
case  when  the  ill-conditioning  of  the  system  occurs  is  if  the  order  of  the  models  become  to  high. 
In  this  case  the  scheme  of  parameter  space  division  can  be  applied  to  construct  the  model,  as 
described  in  section  1 .4.4. 

Example.  To  show  the  benefits  of  adaptive  sampling  over  a  uniform  rectangular  grid  a  com¬ 
parison  of  mentioned  techniques  is  presented.  The  model  order  was  constant  and  selected  as 
Vs\  =  [2  2  2  2].  Model  accuracy  for  different  density  D  of  grid  is  presented  in  table  1.5,  figure 

1.5  shows  the  histograms  of  the  model  error.  To  compute  the  histograms  a  iris  response  on 
10000  random  data  points  was  computed  in  electromagnetic  solver  and  compared  with  model 
response.  The  error  was  computed  as  direct  difference  between  model  and  simulation  response 
A  =  ||S(X)  —  S(X_)  || .  Figure  1.5  also  presents  a  cumulative  distribution  function  of  the  error 
computed  on  the  same  10000  points,  that  shows  which  percent  of  samples  has  accuracy  equal 
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Figure  1.5:  Histogram  and  cumulative  distribution  function  of  model  error  for  different  grid 
density:  a)  D  =  4,  b)  D  —  5,  c)  D  —  6,  d)  D  =  7 


Table  1.5:  Model  accuracy  vs.  density  of  rectangular  grid 


D 

M 

A max  [dB] 

A mean  [dB] 

4 

256 

-26.54 

-40.66 

5 

625 

-27.67 

-42.71 

6 

1296 

-27.97 

-43.46 

7 

2401 

-28.03 

-43.56 

or  better  than  given  error. 

As  expected,  the  more  dense  is  the  rectangular  grid,  the  better  accuracy  of  the  computed 
model  can  be  observed.  However,  comparing  the  improvement  of  the  model  accuracy  to  the 
number  of  support  points,  it  have  to  be  noted  that  increase  of  the  grid’s  density  is  not  an  efficient 
way  to  improve  the  model  quality. 

To  compare,  the  same  model  was  created  with  application  of  adaptive  sampling  technique. 
An  initial  sparse  grid  resolution  was  D  =  3  (the  grid  had  only  81  points).  The  initial  data  set 
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Figure  1.6:  Histogram  and  cumulative  distribution  function  of  model  error  with  adaptive  sam¬ 
pling  technique 


Table  1.6:  Model  accuracy  using  adaptive  sampling  technique 


M 

A max  [dB] 

A mean  [dB] 

100 

-28.15 

-37.85 

125 

-29.55 

-39.41 

150 

-28.49 

-39.91 

175 

-29.95 

-40.44 

250 

-29.51 

-40.50 

was  appended  with  application  of  adaptive  sampling.  The  model  accuracy  vs.  number  of  points 
added  is  depicted  in  table  1.6.  It  can  be  observed  that  a  addition  of  data  set  of  19  points  gives 
a  smaller  maximum  model  error  than  a  rectangular  grid  consisting  of  2401  points,  and  further 
application  of  adaptive  sampling  improves  the  model  accuracy.  The  adaptive  model  constructed 
on  175  data  points  has  a  similar  mean  error  as  the  model  created  on  a  rectangular  grid  of  256 
data  points  and  maximum  error  is  reduced  by  half.  However,  if  the  adaptive  sampling  proce¬ 
dure  continues  with  the  same  model  orders,  supplementing  of  the  data  set  improves  the  model 
accuracy  marginally. 

QR-update.  The  described  procedure  requires  one  to  update  both  models  each  time  a  support 
point  is  added.  It  means,  that  one  have  to  recompute  the  TLS  solution  of  interpolation  problem 
in  every  iteration.  To  reduce  the  numeric  cost  of  such  operation  an  (TK-update  procedure  can  be 
used  [19].  Assume  that  matrix  C  has  a  factorization  C  =  Q-  R,  where  Q  is  orthogonal  and  R  is 
upper  triangular.  Addition  of  a  single  support  point  appends  a  vector  wT  to  the  matrix  C.  As  an 
result  one  obtains  updated  matrix: 


(1.9) 
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Additionally,  one  can  notice  that: 


diag(l,QT)  ■  C  — 


=  H 


(1.10) 


where  H  is  an  upper  Hessenberg  matrix.  It  is  possible  to  apply  a  set  of  n  subsequent  Givens 
rotations  that  transform  H  to  upper  triangular  form: 

Rl=JTnJTn-\  (l.H) 


Once  the  Givens  rotations  are  known,  the  matrix  Q\  can  be  computed  as: 

Qi=diag(l,Q)  JiJ2  ...  Jn  (1.12) 

Matrices  R\  and  Q\  form  a  QR  factorization  of  matrix  C  =  Q\  ■  R\. 

The  full  QR  factorization  form  scratch  is  algorithm  of  complexity  0(N3),  while  update  of 
the  existing  Q  and  R  matrices  is  0(N 2)  algorithm. 


1.4.3  Order  selection 

The  application  of  adaptive  sampling  using  constant  model  orders  is  not  sufficient  to  automated 
model  construction.  The  key  element  in  investigated  interpolation  procedure  is  a  development 
of  efficient  technique  of  model  order  selection. 

The  adaptive  sampling  procedure  uses  two  different  models  which  are  iteratively  compared 
to  each  other.  The  models  S \  and  S2  are  described  with  two  vectors  V51  and  Vs2.  To  ensure  the 
best  performance  of  the  technique  the  initial  models  have  low  complexity  and  an  algorithm  of 
automated  order  selection  is  applied.  The  minimal  model  order  considered  is  V5i(l  :  N)  =  2 
and  is  assumed  that  model  S2  has  a  higher  order  allowed  than  model  5] . 

Orders  of  initial  models 

The  initial  models  should  be  of  a  similar  order.  In  fact,  the  more  both  models  differ  each  other, 
the  more  data  points  is  needed  the  algorithm  to  converge.  However  it  marginally  influence  the 
accuracy  of  final  models,  because  the  lower  order  model  limits  the  accuracy. 

To  prove  this  a  models  of  the  same  rectangular  iris  as  described  before  were  constructed 
applying  adaptive  sampling  with  different  pairs  of  models  order  selected: 

•  VSi  =  [2  2  2  2],  V52  -  [3  2  2  2] 

•  VSi  =  [2  2  2  2],  V52  =  [3  3  22] 

•  Vsi  =  [2  2  2  2],  V52  =  [3  3  3  2] 

The  initial  data  set  was  appended  with  100  points  by  adaptive  sampling.  Figure  1.7  shows  the 
histogram  and  cumulative  distribution  function  of  error  of  models  low  and  high  order.  It  can  be 
seen  that  the  accuracy  is  mainly  function  of  model  order.  In  practical  computations  the  more 
both  models  differ  each  other,  the  more  support  points  is  needed  to  reach  the  required  accuracy 
£0.  To  obtain  the  best  performance  (accuracy  comparing  to  number  of  support  points)  the  small¬ 
est  difference  between  both  models  is  required.  We  propose  to  use  V51  (1  :  N)  —  [2  2  ...  2]  and 
1/52(1  :  N)  —  [3  2  ...  2]  and  improve  the  accuracy  with  efficient  adaptive  model  order  selection. 
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Figure  1.7:  Histograms  and  cumulative  distribution  function  of  model  error  for  different  pairs 
of  models  used  in  the  adaptive  sampling 


Adaptive  model  order  selection 

The  whole  procedure  starts  with  a  sparse  rectangular  grid  of  support  points  and  low  order  mod¬ 
els.  The  technique  monitors  the  error  between  two  models  created  in  an  iterative  way  with 


cumulative  distribution  function  cumulative  distribution  function  cumulative  distribution  function 
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Table  1.7:  Model  accuracy  and  number  of  support  points  vs.  difference  between  models. 


Model 

A max  [dB] 

A mean  [dB] 

Vsi  -  [2  2  2  2] 
V52  =  [3  22  2] 

-31.32 

-34.37 

-40.70 

-48.15 

Vsi  =  [222  2] 
V52  =  [3  3  2  2] 

-31.90 

-35.56 

-41.12 

-47.54 

Vsi  =  [222  2] 
V52  =  [3  3  3  2] 

-30.91 

-34.60 

-40.76 

-47.58 

adaptive  sampling.  The  behavior  of  the  error  is  a  basic  indicator  if  the  model  order  should  be 
increased.  It  was  presented  previously  that  subsequent  addition  of  support  points  improves  the 
model  quality  until  the  stagnation  phase.  The  basic  method  to  detect  this  is  to  observe  the  num¬ 
ber  of  iterations  without  error  improvement.  In  our  tests  we  decided  that  number  from  range 
2n~ 1  up  to  2n  iterations  without  improvement  is  a  good  indicator  that  the  current  order  (de¬ 
scribed  by  vector  V )  is  too  low  and  has  to  be  increased.  To  make  the  algorithm  more  efficient 
only  these  elements  of  vector  V  which  change  make  the  biggest  reduction  of  error  £  are  raised. 

To  show  the  efficiency  of  the  technique  the  waveguide  iris  model  was  created.  The  initial 
grid  with  D  —  3  was  computed  and  subsequent  data  points  were  added  with  adaptive  sampling. 
The  set  of  59  points  was  added  until  the  stagnation  was  detected.  The  model  order  was  increased 
and  the  next  9  points  was  added  until  the  difference  between  model  decreased  below  the  required 
level  0.003.  The  mean  error  of  created  model  is  -49.5dB  and  maximum  error  is  -33.02dB.  The 
histogram  of  model  error  and  cumulative  distribution  function  are  depicted  on  figure  1.8.  It  can 
be  seen  that  proposed  technique  allowed  to  generate  the  model  in  complete  automated  way  with 
good  accuracy. 


~o 

CD 


Figure  1.8:  Histogram  and  cumulative  distribution  function  of  model  error  with  adaptive  sam¬ 
pling  technique  and  automated  order  selection 
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In  some  cases  such  procedure  is  not  sufficient  to  construct  a  single  model  of  device’s  re¬ 
sponse  due  the  ill-conditioning  of  the  interpolation  problem.  It  may  occur  if  the  model  orders 
are  too  high  and  in  this  case  the  technique  of  parameter  space  division  is  applied. 

1.4.4  Division  of  parameter  space 

Division  of  parameter  space  is  essential  if  the  response  of  complex  device  is  modelled  (for 
example  has  several  resonances)  and/or  superb  accuracy  is  requested.  In  such  cases  it  might 
be  impossible  to  construct  a  single  rational  model  which  covers  whole  parameter  space  and 
assures  desirable  accuracy.  To  overcome  the  problem  an  automated  technique  of  parameter 
space  division  was  developed. 

The  most  important  issue  is  to  develop  criteria  of  space  division.  In  proposed  technique  two 
cases  can  cause  the  division  of  parameter  space  during  the  adaptive  sampling  procedure: 

•  The  size  of  the  problem  becomes  too  big  to  be  efficiently  solved 

•  Further  increasing  of  model  orders  leads  to  ill-conditioned  interpolation  problem 

If  one  of  the  conditions  if  fulfilled,  the  adaptive  sampling  stops  and  the  variance  analysis 
for  each  the  model  parameters  is  performed.  The  smaller  is  variance  of  distribution  of  samples 
connected  with  parameter  v,  the  more  data  points  are  concentrated  around  mean  value  of  xt. 
High  concentration  of  data  points  in  some  area  suggests  that  in  that  place/dimension  the  model 
is  poor,  therefore  that  dimension  is  selected  to  be  divided.  Algorithm  creates  two  smaller  sub¬ 
spaces  with  division  of  range  of  selected  parameter  into  halves.  Such  procedure  is  implemented 
as  recursive  one. 

To  illustrate  the  proposed  algorithm  a  simple  case  of  two-variate  function  S(jci,jc2)  is  pre¬ 
sented  on  figure  1.9.  Figure  shows  that  initial  parameter  space  was  divided  into  three  non- 
overlaping  smaller  subspaces.  The  generalization  to  N-dimensions  is  straightforward. 

To  show  the  robustness  of  the  technique,  it  was  used  to  create  very  accurate  model  of  waveg¬ 
uide  iris.  The  required  accuracy  of  model  was  established  as  £o  =  0.001  (-60dB).  The  procedure 
started  from  sparse  grid  of  81  support  points  and  adaptive  sampling  and  order  selection  reduced 
the  error  level  to  value  0.002.  Further  increasing  of  models  order  resulted  in  ill-conditioning  of 
the  problem.  It  is  worth  to  notice,  that  it  would  be  the  maximum  accuracy  possible  to  obtain 
without  application  of  parameter  space  division. 

To  meet  requested  accuracy  the  initial  parameter  space  shown  in  table  1.1  was  sequentially 
divided  into  three  subspaces,  presented  in  table  1.8.  At  first  the  range  of  width  of  iris  was 
divided,  then  in  one  of  the  subspaces  the  frequency  range  was  divided.  For  each  subspace  an 
independent  model  of  device  response  was  created.  The  histogram  and  cumulative  distribution 
function  of  model  error  are  presented  on  figure  1.10.  90%  of  samples  have  accuracy  better  than 
-50dB.  The  mean  error  of  model  is  -55.97dB  and  maximum  error  drops  to  -38.55dB.  The  total 
number  of  support  points  is  M=460. 

1.4.5  Merging  submodels 

The  main  disadvantage  of  space  division  is  a  non-smooth  response  of  the  models  at  point  of 
connection  of  their  domains.  Since  it  is  impossible  to  impose  the  continuity  conditions  directly 
into  model  computation  algorithm,  this  problem  has  to  be  solved  separately  during  computation 
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subspace  I 
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(model  not  completed) 


space  division 


subspace  II 
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A1 

subspace  III 
(model  completed) 


A1 

space  division 

Figure  1.9:  Sample  space  division  in  two-dimensional  case 


of  the  model  response.  The  problem  is  illustrated  in  figure  1.11,  that  shows  a  plot  of  Su  param¬ 
eter  computed  as  an  response  of  two  models  that  cover  this  frequency  range.  The  discontinuity 
of  the  response  can  be  observed  directly  at  the  point  where  the  parameter  space  was  divided. 

In  most  of  model  applications  the  presence  of  response  discontinuity  is  not  an  important 
issue.  It  is  possible  to  perform  a  successful  design  using  such  non-smooth  model  even  when 


Table  1.8:  Parameter  ranges  for  three  subspaces  created  with  automated  parameter  division 
scheme 


Parameter 

Model  I 

Model  II 

Model  III 

frequency 

width 

height 

thickness 

11.855GHz-  18.02GHz 
6.32mm-  11.06mm 
4.74mm  -  7.899mm 
0.2mm  -  2mm 

11.855GHz-  14.9375GHz 

1 1.06mm  -  15.8mm 
4.74mm  -  7.899mm 
0.2mm  -  2mm 

14.9375GHz-  1 8.02GHz 

1 1.06mm  -  15.8mm 
4.74mm  -  7.899mm 
0.2mm  -  2mm 

CHAPTER  1.  SURROGATE  MODELS 


18 


~o 

CD 


Figure  1.10:  Histogram  and  cumulative  distribution  function  of  model  error  in  case  of  auto¬ 
mated  space  division  used 


the  optimization  of  the  structure  is  involved.  However,  if  one  needs  a  smooth  response  in 
whole  parameter  space,  it  is  possible  to  compensate  the  discontinuity.  One  can  use  a  cubic 
spline  interpolation  procedures  in  the  area  of  model  connection,  as  presented  in  figure  1.11. 
In  this  approach,  a  model  response  in  area  of  the  model  connection  is  computed  from  cubic 
spline  interpolation,  generated  from  six  points  located  near  to  the  model  connection  (3  points 
from  each  model  are  taken  into  account).  Application  of  cubic  splines  gives  as  result  smooth 
response  with  continuous  first  derivative  of  the  response.  Additionally  it  is  fast  and  easy  to 
implement. 

1.4.6  Models  of  multi-port  components 

The  technique  described  in  the  previous  paragraphs  can  be  used  to  construct  a  model  of  single 
scattering  parameter  versus  frequency  and  structure  dimensions.  In  practice,  an  engineer  uses 
multi-port  components,  which  are  described  with  scattering  matrix  that  contains  several  scatter¬ 
ing  parameters.  To  create  a  complete  model  of  such  a  device  all  the  elements  of  the  scattering 
matrix  should  be  modelled  in  an  independent  way.  However,  to  speed  up  this  process  signif¬ 
icantly,  the  successive  model  can  utilize  the  results  of  electromagnetic  simulations  that  were 
already  performed.  In  such  a  case,  at  the  beginning  of  model  creation  the  sparse  grid  is  used 
to  create  the  model  of  first  scattering  parameter  with  adaptive  sampling  and  order  selection.  At 
this  stage  the  results  of  simulations  of  all  scattering  parameters  are  stored.  Once  the  procedure 
converges,  all  the  stored  data  points  can  be  used  to  start  the  generation  of  model  subsequent 
scattering  parameter.  Each  time  the  modelling  of  subsequent  scattering  parameter  is  started,  a 
test  for  initial  model  orders  can  be  performed.  The  test  generates  several  models  with  increasing 
orders  and  evaluates  the  biggest  mismatch  between  chosen  pairs.  The  orders  of  a  model  pair 
with  the  smallest  mismatch  are  used  as  initial  orders  for  adaptive  model  construction  scheme. 

The  presented  formulation  of  multi-port  device  model  construction  can  also  be  used  to  create 
models  of  multi-mode  devices,  as  discussed  in  [2]. 
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1.4.7  Complete  algorithm  -  flow  chart 

The  flow  chart  of  proposed  algorithm  is  presented  on  figure  1.12.  In  the  main  loop  an  adaptive 
sampling  of  the  parameter  space  is  performed.  In  this  loop  the  condition  for  increasing  of  model 
orders  is  checked  and  detection  of  ill-conditioning  is  performed. 


Figure  1.11:  Results  of  proposed  technique  of  model’s  response  discontinuity  compensation: 
—  Non-smooth  model  response,  —  spline  compensated  model  response,  points  used  for 
spline  computation 
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Figure  1.12:  Flow  chart  of  complete  algorithm 
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1.5  State-of-the-art  examples 

The  illustrate  the  flexibility  of  proposed  technique  models  of  two  very  complex  planar  devices 
were  created.  In  both  cases  the  structures  were  simulated  using  a  commercial  electromagnetic 
solver  Momentum.  These  examples,  along  with  the  waveguide  iris  model  presented  in  sec¬ 
tion  1.4.4,  shows  that  with  proposed  technique  it  is  possible  to  create  a  surrogates  of  complex 
microwave  devices. 

1.5.1  Spiral  inductor  in  SiGe  BiCMOS  technology 

Structure  overview  of  an  octagonal  spiral  inductor  is  presented  on  figure  1.14.  Figure  1.15 
shows  a  three  dimensional  view  of  the  structure  along  with  current  visualization  on  the  surface 
of  the  inductor.  The  inductor  consists  of  2.5  loop  spiral  with  an  uniform  strip  at  the  top  layer 
and  metal  bridge  at  lower  layer.  The  structure  is  described  with  three  geometric  parameters: 
strip  width  w,  gap  width  s  and  inner  spiral  radius  R.  The  model  of  such  structure  was  computed 
in  frequency  range  from  DC  to  10GHz.  Three  scattering  parameters  were  modelled:  sn,  521 
and  S22 ■  A  range  of  model  parameters  is  presented  in  table  1.9.  Selected  parameter  range 
corresponds  to  approximate  inductor  die  area  range  from  150 fjm  x  150 fjm  up  to  470/rm  x  470/rm. 
The  required  tolerance  was  set  as  £=0.001,  and  the  procedure  needed  285  support  points  to  build 
the  models.  The  time  of  single  analysis  was  from  3  to  5  minutes  on  a  1.5GHz  PC  and  depends 
of  structure  size. 

Then  a  set  of  500  randomly  distributed  data  points  in  the  model  domain  was  generated.  The 
set  was  used  to  verify  the  accuracy  of  the  model.  The  mean  error  computed  over  the  set  was 
— 67. 3dB,  while  maximum  error  reaches  — 55dB.  Figure  1.16  shows  a  histogram  and  cumulative 
distribution  function  of  the  model  error  of  5n.  It  can  be  seen,  that  90%  of  samples  has  an  error 
below  — 63dB.  Additionally  the  histogram  and  cumulative  distribution  function  of  error  of  Q 
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Figure  1.13:  Silicon  substrate  as  simulated  in  ADS  Momentum 
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w  s 


Figure  1.14:  Structure  and  dimensions  of 
modelled  octagonal  inductor. 


Figure  1.15:  Three  dimensional  field  visu¬ 
alization  of  modelled  inductor. 


factor  computed  from  model  response  is  presented  in  Fig.  1.1 7.  Despite  good  accuracy  of  the 
models  of  scattering  parameters,  the  error  of  extracted  Q  factor  can  be  high.  However  the  high 
error  appears  very  near  of  the  inductor  resonance  (scattering  matrix  is  close  to  unitary  matrix) 
and  in  this  point  the  parasitic  capacitance  of  the  inductor  dominates.  In  fact  this  point  of  inductor 
operation  is  out  of  applications. 


Table  1.9:  Range  of  input  parameters  of  spiral  inductor  model 


Parameter 

Range 

frequency  (/) 
strip  width  (w) 
gap  width  (5) 
inner  radius  ( R ) 

0GHz  -  10GHz 
10 fim  -  25 /urn 

5 /urn  -  20 /urn 

30 /um  -  100 /jtn 

1.5.2  Interdigitated  capacitor  in  MCM-D  technology 

A  next  example  is  an  interdigital  capacitor  made  in  MCM-D  technology.  Thin  film  MCM-D 
technology  has  many  advantages  over  the  traditional  hybrid  technologies.  It  assures  high  preci¬ 
sion  components,  repeatability  of  manufacturing  complex  microwave  structures  and  integration 
of  analog  and  digital  circuits. 

A  layout  of  an  interdigitated  capacitor  made  in  MCM-D  technology  is  shown  on  figure  1.18. 
A  model  of  six  variables  (frequency  and  five  geometric  dimensions)  was  created.  The  ranges 
of  the  model  parameters  are  presented  in  table  1.10.  The  substrate  parameters  are  as  follows: 
thickness  45/um,  dielectric  permittivity  er  —  2.65,  dielectric  losses  tan 8  =  0.002,  metal  thickness 
3/um,  metal  conductivity  8  =  4.525  ■  107^  (gold). 

The  requested  modelling  error  was  £  =0.001  =-60dB.  The  procedure  started  using  a  sparse 
grid  with  D= 3  (729  support  points)  and  initial  orders  V  —  [2  2  2  2  2  2].  Adding  of  1402  data 
points  with  increasing  of  orders  to  V  —  [3  4  3  4  4  4]  gives  the  mismatch  between  both  models 
below  the  requested  0.001 .  To  verify  the  model  accuracy  a  set  of  3000  randomly  distributed  data 
points  were  computed  using  an  electromagnetic  simulator  and  compared  with  model  response. 
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Figure  1.16:  Histogram  and  cumulative  distribution  function  of  Sn  model  of  octago¬ 
nal  inductor. 


Figure  1.17:  Histogram  and  cumulative  distribution  function  of  (7- factor  error  of  ra¬ 
tional  model  of  octagonal  inductor. 


CHAPTER  1.  SURROGATE  MODELS 


24 


Table  1.10:  Range  of  input  parameters  of  interdigitated  capacitor  model 


Parameter 

Range 

frequency  (/) 
input  line  width  (wi) 
cap.  line  width  (W2) 
finger  length  (L) 
finger  line  width  (w3) 
gap  width  (g) 

10GHz  -  90GHz 
50 jum  -  150 jum 
lOjum  -  25 jim 
100 nm  -  250 /um 
10 /an  -  20 fim 

10 /an  -  25 fan 

The  histogram  and  cumulative  distribution  function  of  model  error  are  depicted  on  figure  1.19. 
The  maximum  error  of  created  model  is  -32.5dB  and  mean  error  reaches  -56.5dB.  In  case  of 
90%  of  test  samples  the  error  is  below  -50dB  and  for  54%  of  samples  the  error  is  below  -60dB. 

1.5.3  Summary 

In  table  1.11  presented  is  a  detailed  comparison  of  described  models.  It  can  be  seen,  that  the 
model  accuracy  strongly  depends  on  complexity  of  the  structure  and  number  of  variables. 


Table  1.11:  Comparison  of  created  models 


Structure 

N 

M 

Space  divisions 

A max  [0B] 

A mean  [0B] 

A90  [dB] 

Octagonal  inductor 

4 

285 

0 

-55 

-67.3 

-63 

Waveguide  iris 

4 

460 

2 

-38.55 

-55.97 

-51 

MCM-D  Capacitor 

6 

2131 

0 

-32.5 

-56.5 

-50 

1.6  Applications 

The  proposed  technique  is  versatile  and  can  be  applied  for  construction  of  surrogate  models  of 
different  kind  of  microwave  structures.  In  fact,  every  structure  that  is  described  with  scattering 
parameters  and  has  a  smooth  response  can  be  modelled.  Therefore  one  can  create  models  of 
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Figure  1.19:  Histogram  and  cumulative  distribution  function  of  model  error  for  interdigitated 
capacitor 


elements  made  in  various  technologies:  waveguide,  microstrip,  coplanar  or  multilayer.  To  show 
the  benefits  of  models  application  two  advanced  examples  are  presented  below. 


1.6.1  Inductor  optimization 


High  Q  inductors  is  one  of  the  problems  in  on-chip  solutions  for  RF  and  microwaves.lt  is  a 
common  problem  of  designers  to  find  an  inductor  that  at  frequency  /  realizes  an  inductance  Lq 
and  in  the  same  time  minimize  the  inductor’s  area  S  and/or  maximize  the  Q.  It  is  an  important 
issue,  since  there  exist  several  designs  of  inductor  with  the  same  value  of  inductance  Lq  and 
different  values  of  quality  factor  and  size.  Several  papers  have  investigated  the  inductor  opti¬ 
mization  problem  [20,  21,  22,  23].  Authors  propose  to  create  a  simple  closed-form  model  of 
inductor  parameters,  like  in  [20,  23],  but  such  models  are  difficult  to  find  for  multivariate  case. 

This  problem  can  be  solved  using  optimization  procedure.  Frequency  dependent  parameters 
of  the  inductor  can  be  evaluated  from  its  admittance  parameters  as  [20] : 


1 

2nf  ■  Im(Y\  i  (/)) 


(1.13) 


MYn(f)) 


(1.14) 


Both  parameters  L  and  Q  of  inductor  depend  of  frequency  /  and  dimensions  of  the  structure 
X  =  [x\  X2  ...  v,v] .  The  procedure  of  inductor  search  can  be  organized  as  a  simple  optimization 
routine  with  goal  function  G  defined  as: 


N 


G(fo,X)  =  \Wo,X)  —T0||  —  Q(fo,X_)  + 

1=1 


(1.15) 
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which  has  to  be  minimized  over  space  of  inductor  geometric  parameters.  The  minimum  of 
the  goal  function  gives  the  desired  value  of  L  for  the  smallest  inductor  and  highest  Q  factor. 
Obviously  such  search  is  time  consuming  if  electromagnetic  simulator  is  involved  to  compute 
the  inductor  parameters.  However,  such  search  can  be  performed  very  fast  using  a  multivariate 
surrogate  model  of  the  inductor.  To  prove  this,  a  model  of  octagonal  inductor  described  in 
section  1.5.1  was  investigated  and  several  optimizations  were  performed  using  the  above  goal 
function.  A  few  examples  of  designed  inductors  are  shown  in  table  1.12.  The  time  of  design  in 
all  cases  is  below  1  second.  It  is  extremely  fast  comparing  to  electromagnetic  approach  since 
the  EM-simulation  of  inductor  response  on  single  frequency  takes  a  few  minutes.  It  has  to  be 
noted  that  the  model  parameter  space  is  continuous  but  a  foundry  production  process  assumes 
limited  resolution.  Therefore  final  dimensions  have  to  be  rounded  to  obey  the  process  design 
rules. 

Table  1.12:  Results  for  an  octagonal  inductor  design:  requested  parameters,  obtained  (2- factor, 
approximated  inductor  size  and  time  of  optimization  (MATLAB) 


Requested  parameters 

Q 

Time  [s] 

Aprox.  size 

fo  =  0.92GHz,  L0  =  2.2nH 

10.4 

0.52 

375/r  m  x  375/j  m 

f0=  1.8GHz,  Lq  =  1.4nH 

12.9 

0.49 

208 /u  m  x  208 /u  m 

fo  =  2.45GHz,  Lq  =  2nH 

16.8 

0.7 

236 /j  m  x  236 /j  m 

f0  =  5.1GHz,  L0=  l.lnH 

17.7 

0.45 

172 /r  m  x  172/j  m 

1.6.2  Waveguide  filter 

Surrogate  models  can  be  applied  to  design  complex  microwave  components.  For  example,  a 
5-th  order  microwave  filter  with  two  dispersive  stubs  was  divided  into  a  separate  discontinuities 
that  have  been  modelled  using  the  proposed  approach,  as  presented  in  figure  1.20.  In  the  same 
figure  a  comparison  of  the  model  response  and  the  electromagnetic  simulator  is  presented.  In 
can  be  seen  that  model  response  is  very  accurate.  What  is  important,  it  takes  about  1  second 
to  compute  model  response,  which  is  much  significantly  faster  comparing  to  over  2  minutes  in 
case  of  electromagnetic  analysis  (mode-matching). 


1.7  Problems  with  grid-based  solvers 

The  most  important  issue  for  a  successful  creation  of  rational  model  based  on  the  electromag¬ 
netic  simulation  is  a  smooth  change  of  simulator  response  due  to  changes  of  structure  geometry. 
It  is  a  common  feature  of  mode-matching  based  simulators,  but  may  cause  problems  if  a  grid- 
based  solvers  (like  MoM  or  FDTD)  are  used.  The  problem  is  illustrated  in  figure  1 .21,  where  the 
Si  i  response  of  a  microstrip  stub  on  the  MCM-D  substrate  versus  width  of  the  stub  is  presented. 
The  structure  is  simulated  at  the  frequency  1GHz  and  re-meshed  each  time  the  structure  dimen¬ 
sions  changes.  Mesh  frequency  is  also  set  as  1GHz.  It  can  be  seen  that  the  transfer  function  is 
not  smooth,  which  indicates  non-physical  response.  A  discontinuity  of  the  response  causes  a 
huge  problem  for  interpolation  of  data  using  rational  functions. 
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f  [GHz] 

Figure  1 .20:  Waveguide  filter  structure  and  comparison  of  model  response  with  results  of  elec¬ 
tromagnetic  simulator  (mode-matching). 


A  way  to  circumvent  this  problem  is  to  simulate  the  structure  with  a  grid  that  is  denser  than 
the  one  resulting  by  considering  the  frequency  alone.  In  the  same  figure  the  response  of  the 
same  device  is  shown  with  mesh  frequency  set  as  100GHz.  In  this  case  the  response  is  smooth 
with  change  of  the  geometry.  A  drawback  of  such  solution  is  an  increased  simulation  time,  due 
to  denser  mesh. 
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Figure  1.21:  Response  of  the  stub  on  the  MCM-D  substrate  in  function  of  stub  width:  —  mesh 
computed  at  1GHz,  —  mesh  computed  at  100GHz 


Chapter  2 

Integration  with  circuit  simulators 


The  second  goal  of  the  grant  was  to  integrate  the  developed  techniques  with  the  industry  stan¬ 
dard  circuit  simulators.  Two  issues  are  of  importance: 

•  Generation  of  support  points 

•  Integration  of  surrogate  models  with  circuit  simulators 


2.1  Obtaining  data  from  planar  simulators 

The  algorithm  of  automatic  model  creation  depends  on  the  results  of  electromagnetic  simula¬ 
tions.  Because  of  continuity  of  parameter  space  and  numerous  support  points,  the  algorithm 
must  be  able  to  create  necessary  structure,  analyze  it  in  a  simulator  and  get  the  results.  Most 
popular  planar  simulators  can  be  used  as  a  source  of  data  samples,  e.g.  Sonnet,  AWR  Mi¬ 
crowave  Office  or  Agilent  Momentum.  Details  of  implementations  will  be  explained  taking  as 
an  example  simple  microstrip  tee  structure,  simulated  on  single  frequency  point  2 GHz,  shown 
on  Fig.2.1.  The  substrate  parameters  are  as  follows:  height  0.25 mm,  dielectric  constant  9.6, 
conductor  conductivity  1  ■  1050S  without  dielectric  losses  and  enclosure  without  cover. 

2.1.1  Sonnet  simulator 

Sonnet  simulator  can  be  run  in  a  batch  mode  from  command  line  taking  as  an  input  a  text  file 
containing  all  information  needed  for  analysis.  This  is  done  by  running  a  command: 

em.exe  input_f ile . son 

Results  of  simulation,  i.e.  scattering  parameters,  are  stored  in  file  log-response.log  in  son- 
data  subfolder.  Simple  text  file  parsing  is  sufficient  to  extract  numeric  data. 

Because  of  the  complexity  of  simulator  and  great  number  of  options,  only  those  relevant 
to  modeling  application  are  described  in  the  next  section.  The  input  file  format  description  is 
based  on  version  10.51  of  Sonnet  released  in  December  2004. 
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Figure  2.1:  An  example  of  microstrip  tee  structure. 


Input  file  format 

The  first  two  lines  of  the  input  file  contain  only  identification  information  irrelevant  to  simula¬ 
tion.  For  example,  following  two  lines  indicate  a  sonnet  project  file  created  in  version  10.51  of 
the  simulator. 

FTYP  S0NPR0J  1  !  Sonnet  Project  File  VER  10.51 

The  rest  of  the  file  is  divided  into  following  sections: 

•  HEADER 

•  DIM 

•  FREQ 

•  CONTROL 


•  GEO 
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Each  section  begins  with  section  jiame  string  and  ends  with  END  section  jiame  string.  For 
example: 

HEADER  . . .  END  HEADER 

HEADER  section 

The  first  section  has  a  similar  purpose  as  the  first  two  lines  of  the  file.  It  contains  information 
not  related  to  the  simulation,  e.g.  information  about  license  and  dates  of  creation,  modification 
of  the  project,  and  can  be  omitted  during  the  input  file  generation  for  modeling.  An  example  of 
this  section  looks  as  follows: 

HEADER 

LIC  SL00000 .101 

DAT  04/11/2006  14:38:49 

BUILT_BY_CREATED  xgeom  10.51  04/11/2006  14:25:39 
BUILT_BY_SAVED  xgeom  10.51 
MDATE  04/11/2006  14:26:38 
HDATE  04/11/2006  14:26:02 
END  HEADER 

DIM  section 

The  DIM  section  must  contain  definitions  of  units  of  several  dimensions.  All  numerical 
values  in  the  input  file  must  correspond  to  the  units  defined.  Definition  of  every  dimension  unit 
takes  one  line  a  has  a  following  format: 

dim_id  dim_unit 

where  dim  jd  is  identification  string  of  the  dimension  and  dim_unit  is  a  string  indicating 
the  unit.  The  table  below  shows  all  the  dimension  which  definitions  of  units  are  needed. 


Table  2.1:  Description  of  dimension  units. 


Dimension 

dimrid 

dim unit 

Frequency 

FREQ 

Hz,  KHz,  MHz,  GHz,  THz,  PHz 

Inductance 

IND 

H,  MH,  UH,  NH 

Fongitude 

FHG 

M,  MM,  UM,  MIF,  IN 

Angle 

ANG 

DEG 

Conductivity 

CON 

/OH 

Capacitance 

CAP 

F,  MF,  UF,  NF,  PF 

Resistance 

RES 

OH 

An  example  of  the  DIM  section  is  shown  below: 

DIM 

FREQ  GHZ 
IND  NH 
LNG  MM 
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ANG 

DEG 

CON 

/OH 

CAP 

PF 

RES 

OH 

END  DIM 

FREQ  section 

Frequency  points  are  defined  in  FREQ  section.  It  is  possible  to  define  them  in  the  same  way 
as  in  GUI  of  the  simulator.  Thus  frequency  point  can  be  defined  as  follows: 

STEP  freq_point 

or 

SIMPLE  low_freq  high_freq  interval 

where  freq_point  is  a  frequency  of  a  single  point,  low_freq  and  high_freq  are  limits  of  a  set 
of  frequency  points  with  interval  between  the  points.  For  example  the  following  line: 

STEP  2.455 

sets  a  frequency  point  at  2.455  units  of  frequency,  and  the  next  line: 

SIMPLE  2.50  2.80  0.1 

defines  31  frequency  points  from  2.50  to  2.80  units  with  an  interval  of  0.1  units.  Other  types  of 
frequency  points  definition  are  irrelevant  for  modeling  purposes. 

CONTROL  section 

This  section  includes  simulation  options.  The  most  important  ones  are: 

SPEED  value 

where  value  indicates  preset  setting  of  mesh  density  influencing  accuracy  and  speed  of  analysis. 
The  value  can  be  0,  1  or  2.  The  lower  the  value  is,  the  more  accurate  and  slower  simulation  is. 
Another  useful  command  in  CONTROL  section  is: 

OPTIONS  -d 

meaning  that  ports  will  be  de-embedded  from  results. 

GEO  section 

This  section  of  the  input  file  contains  all  information  about  geometry  of  the  structure,  in¬ 
cluding  box,  dielectric  and  metallization  layers  parameters.  First,  the  top  and  the  bottom  metals 
must  be  defined.  It  can  be  chosen  from  three  predefined  types: 

•  Lossless,  models  a  perfect  conductor 
TMET  "Lossless"  0  SUP  0  0  0  0 

•  Load,  models  a  perfect  matched  waveguide  load 


TMET  "WG  Load"  0  WGLOAD 
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•  Free  Space,  which  removes  the  top  or  bottom  cover 

TMET  "Free  Space"  0  FREESPACE  376.7303136  000 

It  is  also  possible  to  set  a  user-defined  metal  as  the  top  or  the  bottom  metal.  For  the  top  metal¬ 
lization  it  can  be  done  placing  following  line: 

TMET  "metal_name"  metal_.no  NOR  metal_cond  metal_cr  metal_thick 

where  metal_name  is  a  unique  name  of  the  metal,  metal_no  is  a  successive  number  of  the 
metal,  i.e.  1  for  the  first  user-defined  metal,  metaLcond  is  the  metal  conductivity,  metaLcr  is 
the  current  ratio  and  metaLthick  is  the  thickness  of  the  metal. 

It  is  required  that  user-defined  metal  must  be  declared  in  a  separate  line  in  the  same  way  as 
TMET  but  with  MET  command,  for  example: 

MET  "Metal2 "  1  NOR  1000  0  0.01 

If  a  circuit  is  symmetric  about  the  center  line  parallel  to  the  X  axis,  then  the  symmetry  can 
be  activated  by  placing  the  following  line: 

SYM 

which  results  in  faster  analysis. 

Specification  of  the  box  containing  the  structure  consists  of  a  line  defining  box  dimensions 
and  lines  containing  dielectric  layers  information.  The  box  dimensions  are  defined  as: 

BOX  met_lay_no  box_x  box_y  cells_x  cells_y  20  0 

where  met  Jay_no  is  a  number  of  metallization  layers,  box_x  and  box_y  are  dimensions  of  the 
box,  cells_x  and  cells_y  are  twice  the  numbers  of  cells  horizontally  and  vertically  respectively. 
For  example  the  following  line: 

BOX  4  320  320  256  256  20  0 

defines  a  box  with  four  layers  of  metallization  (five  dielectric  layers)  with  dimensions  320x320 
units  divided  on  128  cells  horizontally  and  vertically. 

Each  dielectric  layer  must  be  defined  in  a  separate  line  of  the  form: 

d_thick  d_const  m_perm  d_loss  m_loss  d_cond  0  "d_name" 

where  cLthick  is  the  thickness  of  the  layer,  dLconst  is  the  dielectric  constant,  m_perm  is  the 
relative  magnetic  permeability,  d  Joss  is  the  dielectric  loss  tangent,  m  Joss  is  the  magnetic  loss 
tangent,  d_cond  is  the  dielectric  conductivity  and  d_name  is  a  unique  name  of  the  layer.  For 
example  the  following  line: 

25  12.9  1  0  0  0  2  "GaAs" 

defines  a  lossless  layer  of  gallium  arsenide  25  units  thick. 

Definition  of  ports  spans  over  four  lines.  The  first  line  begins  port  declaration: 


PORI  port_type 
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where  port_type  can  be  STD  for  a  standard  port  or  AGND  for  an  autogrounded  port.  The 
second  line  binds  the  port  with  a  specific  polygon  of  metallization: 

POLY  poly_id  1 

where  poly Jd  is  a  unique  number  of  the  polygon  to  which  the  port  is  adjacent..  The  third  line 
defines  an  edge  the  box  the  port  is  adjacent  to.  For  example: 

3 

means  that  the  port  is  adjacent  to  left  edge  of  the  box.  The  other  values  are:  0  for  top  edge,  1 
for  right  edge  and  2  for  bottom  edge.  The  last  line  of  port  declaration  declares  the  number  of 
the  port,  the  port  impedance  and  its  position: 

port_no  port_res  port_react  port_ind  port_cap  port_x  port_y 

where  port_no  is  a  successive  number  of  the  port,  port_res  and  port_react  are  real  and  imagi¬ 
nary  parts  of  the  port  impedance,  port_ind,  port_cap  are  inductance  and  capacitance  of  the  port 
and  port_x  and  port_y  are  coordinates  of  the  port.  For  example  the  following  line: 

2  50  0  0  0  320  160 

defines  a  second  port  with  purely  resistive  impedance  of  50  units  at  position  x=320  and  y=160. 
Reference  planes  can  be  moved  away  from  the  edge  of  the  box  by  placing  a  line: 

DRP 1  rp_pos  FIX  rp_dist 

where  rp_pos  indicates  the  edge  of  the  box  of  reference  plane  displacement  (can  be  LEFT, 
RIGHT,  TOP  or  BOTTOM),  rp_dist  is  a  distance  of  the  reference  plane  from  the  specified 
edge.  For  example,  the  following  line: 

DRP 1  RIGHT  FIX  1.20 

sets  reference  place  for  ports  on  the  right  edge  at  a  distance  of  1.20  units  from  the  edge. 

Shapes  of  metallization  or  dielectric  are  drawn  in  a  form  of  polygons.  Description  of  a  poly¬ 
gons  begins  with  a  line  indicating  the  total  number  of  polygons  in  the  structure.  For  example: 

NUM  2 

indicates  there  are  two  polygons  to  be  defined  below.  Declaration  of  a  single  polygon  begins 
with  the  following  line: 

layer_no  poly_vert  metal_.no  N  poly_id  sub_xmin  sub_ymin  sub_xmax 
sub_ymax  000  egde_mesh 

where  layer_no  is  a  metallization  layer  number,  poly_vert  is  a  number  of  vertices  in  the 
polygon  (including  a  duplicate  vertex  at  the  end  of  the  list),  metal_no  is  metal  number  (equiv¬ 
alent  to  metal  number  in  user-defined  metals),  poly  Jd  is  an  unique  polygon  id,  sub_xmin, 
sub_ymin,  sub_xmax  and  sub_ymax  are  subsectioning  controls,  egdejnesh  turns  on/off  edge 
mesh  and  can  be  either  Y  or  N. 

Dimensions  of  the  polygon  are  specified  by  defining  the  position  of  all  vertexes  in  the  poly¬ 
gon.  The  position  of  each  vertex  is  described  by  the  line: 
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poly_x  poly_y 

where  poly_x  and  poly_y  are  positions  relative  to  a  coordinate  system  of  the  structure.  Note  that 
the  last  vertex  must  be  the  same  as  the  first  one. 

For  example  block  below: 

0  5  -1  N  12  1  1  100  100  0  0  0  Y  320  147.5  320  172.5  220  172.5  220 
147.5  320  147.5  END 

defines  a  rectangle  in  the  first  metallization  layer  with  id  12  with  dimensions  100x25  units  and 
bottom  left  corner  at  the  position  x=220  and  y=147.5  units. 

Example 

The  structure  presented  on  Fig. 2.1  would  be  described  by  a  file  which  most  important  sections 
are: 


•  FREQ  section  for  a  single  frequency  point  at  2 GHz 

FREQ 

SIMPLE  2.0 
END  FREQ 

•  GEO  section  with  box  and  dielectric  layers  parameters 


GEO 

TMET  "Free  Space"  0  FREESPACE  376.7303136  000 
BMET  "Lossless"  0  SUP  0  0  0  0 
BOX  1  7  6.2  140  124  20  0 

30  1  1  0  0  0  0  "Air" 

0.25  9.6  1  0  0  le+050  0  "substrate" 


...  ports  definitions  ... 

PORI  STD 
POLY  1002  1 
3 

1  50  0  0  0  0  3.1 
PORI  STD 

POLY  1003  1 

2 

3  50  0  0  0  3.5  6.2 
PORI  STD 
POLY  1005  1 
1 

2  50  0  0  0  7  3.1 


...  and  metallization  polygons 
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NUM 

4 

0  5 

- 

1  N 

1002 

1 

1 

0  2, 

.5 

2.5 

2 

.5 

2.5 

3 

.7 

0  3, 

.7 

0  2, 

.5 

END 

0  5 

- 

1  N 

1003 

1 

1 

2.5 

3 

.7 

4.5 

3 

.7 

4.5 

6 

.2 

2.5 

6 

.2 

2.5 

3 

.7 

END 

0  5 

- 

1  N 

1005 

1 

1 

4.5 

2 

.8 

7  2, 

.8 

7  3, 

.4 

4.5 

3 

.4 

4.5 

2 

.8 

END 

0  7 

- 

1  N 

1006 

1 

1 

2.5 

3 

.7 

2.5 

2 

.5 

3.5 

2 

.5 

3.5 

2 

.8 

4.5 

2 

.8 

4.5 

3 

.7 

2.5 

3 

.7 

END 

END 

GEO 

100  100  0  0  0  Y 


100  100  0  0  0  Y 


100  100  0  0  0  Y 


100  100  0  0  0  Y 


2.1.2  AWR  Microwave  Office 

Microwave  Office  implements  a  COM-compatible  interface  allowing  easy  communication  with 
other  software.  All  functions  related  to  creating  EM  structures  within  simulator’s  GUI  have 
COM  counterparts.  Thus  creating  the  structure  is  straightforward  and  consists  of  the  following 
steps  presented  in  pseudocode  taking  as  an  example  the  structure  from  Fig.  2.1: 

1 .  Starting  new  project 

proj  =  MW0.New(  'eoard'  ) 


2.  Creating  an  empty  EM  structure 
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Figure  2.2:  View  of  the  example  structure  in  Sonnet  window. 


emstructure  =  proj .EMStructures.Add(  'mstrip  tee'  ) 


3.  Setting  substrate  -  when  creating  the  EM  structure,  there  are  two  default  material  lay¬ 
ers.  Thus  no  new  layers  must  be  created  in  this  case  -  only  the  existing  ones  need  to  be 
modified.  For  setting  the  substrate  parameters: 

layer2  =  MaterialLayers . Item (2) 
layer2 . Thickness  =  0.25 
layer2 .DielectricConstant  =  9.6 
layer2 . LossTangent  =  0.0 


and  for  the  air  layer: 

layerl  =  MaterialLayers . Item(l) 
layerl . Thickness  =4.0 
layerl .DielectricConstant  =  1.0 
layerl . LossTangent  =  0.0 


Setting  the  enclosure  requires  the  following  commands: 

emstructure . Enclosure . XDimension  =  0.0070 
emstructure . Enclosure . YDimension  =  0.0062 
emstructure . Enclosure . XDivisions  =  70 
emstructure . Enclosure . YDivisions  =  62 
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emstructure . Enclosure . Height  =4.0 

emstructure . Enclosure . EnclosureTop  =  mwBMT_ApproximateOpen 


4.  Drawing  metallization  shapes 


emstructure . Shapes . AddRectangle 
emstructure . Shapes .AddPolygon (  [ 


emstructure . Shapes .AddRectangle 
emstructure . Shapes .AddRectangle 


0,  0.0025,  0.0025,  0.0012  ) 
0.0010  0.0025; 

0.0025  0.0020; 

0.0045  0.0025; 

0.0030  0.0045; 

0.0034  0.0040; 

0.0035  0.0034; 

0.0050  0.0035; 

0.0037  0.0060; 

0.0025  0.0037  ]  ) 

0.0025,  0,  0.0020,  0.0025  ) 
0.0045,  0.0028,  0.0025,  0.0006 


5.  Adding  ports 

emstructure 

emstructure 

emstructure 


Shapes . AddPort ( 
Shapes .AddPort ( 
Shapes .AddPort ( 


0.0000,  0.0025, 
0.0070,  0.0028, 
0.0025,  0.0000, 


0.0000,  0.0031,  1  ) 
0.0070,  0.0034,  1  ) 
0.0045,  0.0000,  1  ) 


6.  Setting  frequency  of  simulation 

pro j .Frequencies .Add (  2.000  ) 


7.  Launching  simulation 

pro j . Simulate ( ) 


8.  Obtaining  simulation  results  -  adding  a  graph: 

graph  =  pro j . Graphs .Add (  'mstrip  tee  s-parameters' ,  4  ) 


setting  its  parameters: 

mes  =  graph. Measurements .Add (  'mstrip  tee',  ' Re (S(l,l))'  ) 
mes  =  graph. Measurements .Add (  'mstrip  tee',  'Im(S(l,l))'  ) 
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and  getting  results 

result  =  mes.YValues; 


Detailed  description  of  Microwave  Office  API  can  be  found  in  [24] . 


2.1.3  ADS  Momentum 

Generation  of  a  structure  in  Momentum  is  not  as  straightforward  as  in  Sonnet  or  Microwave 
Office  because  of  numerous  files  generated  by  Momentum  engine.  Even  so  shapes  of  metalliza¬ 
tion,  stored  in  proj.  file,  and  substrate  information,  stored  in  proj.sub,  are  simple  text  files. 


Example 


The  ports  and  metallization  shapes  are  defined  in  proj.  file: 


UNITS  MM, 10000; 


ADD 

Nl 

:  Fl .  0 

:  SO 

:  T1003 

:  D 

'1, 

50, 

0, 

O' 

ADD 

N1 

:F1 . 0 

:  SO 

:  T1003 

:  D 

'  2 , 

50, 

0, 

O' 

ADD 

Nl 

:  Fl .  0 

:  SO 

:  T1003 

:  D 

'  3, 

50, 

0, 

O' 

ADD  PI 


0,  3.1; 
7,  3.1; 
3.5,  0; 


0.000000,3.700000 

2.500000,3.700000 

2.500000,2.500000 
0.000000,2.500000 
0.000000,3.700000; 
ADD  PI 


2.500000.2.500000 

4.500000.2.500000 
4.500000, 0.000000 
2.500000, 0.000000 
2.500000,2.500000; 

ADD  PI 

4.500000.3.400000 
7.000000,3.400000 
7.000000,2.800000 
4.500000,2.800000 
4.500000,3.400000; 

ADD  PI 

3.500000.3.400000 

4.500000.3.400000 
4.500000,2.800000 

4.500000.2.500000 

2.500000.2.500000 

2.500000.3.700000 

3.500000.3.700000 
3.500000,3.400000; 
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The  box  and  dielectric  layers  are  described  in  proj.sub  file: 

UNITS  METRE 

B0TT0MPLANE  IMPEDANCE  0  0 
TOPPLANE  OPEN 

LAYERS 

0  THICKNESS  INFINITY 
PERMITTIVITY  VALUE  1  0 
PERMEABILITY  VALUE  1  0 

r 

1  THICKNESS  0.00025 

PERMITTIVITY  VALUE  9.60 
PERMEABILITY  VALUE  1  0 
STRIP 


Figure  2.3:  View  of  the  example  structure  in  ADS  Momentum  window. 
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2.2  Integration  with  circuit  simulators 

Surrogate  models  created  using  the  presented  technique  can  be  incorporated  into  industry  stan¬ 
dard  circuit  simulators  like  Agilent’s  Advanced  Design  System  or  AWR’s  Microwave  Office. 

2.2.1  ADS  Schematic 

Advanced  Design  System  allows  one  to  create  user-defined  models  using  one  of  the  following 
method: 

1 .  Model  Composer 

2.  Advanced  Model  Composer 

3.  User-Compiled  Model 

Model  Composer  and  Advanced  Model  Composer  create  automatically  models  of  layout  ele¬ 
ments  that  can  be  used  in  linear  simulations.  They  use  polynomial  interpolation  with  reflective 
exploration  technique.  Elements  that  can  be  modeled  using  Model  Composer  are  limited  to  a 
set  of  typical  microstrip  discontinuities,  e.g.  bends,  corners,  crosses  or  gaps  on  a  single  non- 
parameterized  substrate.  The  substrate  can  be  parameterized  in  Advanced  Model  Composer, 
but  for  both  techniques  recommended  number  of  continues  parameters  is  two,  thus  limiting 
modeled  elements  to  very  simple  ones. 

User-Compiled  Model  function  does  not  create  any  models  itself,  but  allows  incorporating 
user-coded  models  of  any  layout  element  and  potentially  removing  the  limitations  of  the  number 
of  parameters.  Linear  user-compiled  models  can  have  up  to  99  external  ports  but  there  is  no 
hard-coded  limitation  regarding  number  of  parameters.  The  making  of  a  model  consists  of 
three  steps: 

1 .  Defining  the  parameters  whose  values  will  be  entered  from  the  schematic 

2.  Defining  the  symbol  and  the  number  of  pins 

3.  Writing  the  C-code 

Definition  of  the  parameters  of  the  model  consists  of  declaring  parameter’s  name,  descrip¬ 
tion,  default  value  and  type.  Validation  of  parameter’s  range  must  be  implemented  by  the  user. 

The  C-code  must  return  Y  matrix  of  the  element  provided  with  element  parameters.  User- 
Compiled  Model  function  automatically  creates  a  template  code  with  a  header  file  for  necessary 
procedures  that  must  be  complemented  by  a  user  code.  Although  the  code  must  conform  to 
ANSI-C  standard,  C++  compiler  is  needed  for  linking  the  entire  program. 

Implementation  details 

Details  of  model  implementation  in  ADS  are  presented  taking  as  an  example  the  MCM-D  ca¬ 
pacitor  model.  The  first  step  is  the  definition  of  model  parameters  depicted  on  Fig. 2. 4.  For 
every  parameter  its  name,  default  value  and  type  are  defined. 

The  next  step  is  to  draw  schematic  symbol  for  the  model  (Fig. 2. 5)  where  the  number  of  pins 
defines  number  of  ports  of  the  model. 
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Figure  2.4:  Parameters  editor  window. 


•S  [  EOARD_PRJ  ]  capacitor  (SYMBOL)  (Schematic):5 


File  Edit  Select  View  Insert  Options  Tools  Layout  Simulate  Window  DynamicLink  DesignGuide  Help 


Figure  2.5:  View  of  schematic  window  when  editing  model  symbol. 


The  main  step  is  writing  the  C-code  for  the  model  calculating  its  Y-parameters.  A  diagram 
of  rational  model  implementation  is  depicted  on  Fig. 2. 6,  where  the  main  part  of  the  implemen¬ 
tation  is  a  function: 

boolean  compute_y (UserlnstDef  *userlnst,  double  omega,  COMPLEX  *yPar) 

The  parameters  typed  in  a  schematic  window  are  passed  into  the  function  in  an  array 


user!nst->pData [] 
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Frequency 

Physical  dimmensions  c=£> 
Substrate  parameters 
etc. 


compute_y() 


Parameters 

i=> 

Rational  function(s) 

S-to-Y  matrix 

validation 

evaluation 

conversion 

Admittance 

parameters 


Figure  2.6:  Diagram  of  model  implementation  in  Agilent’s  ADS. 


and  are  easily  accessible  through  macros  defined  in  the  header  file.  Frequency  of  simulation  is 
passed  in  omega  parameter.  The  yPar  parameter  is  a  pointer  to  an  array  of  complex  admittance 
parameters  that  the  function  must  return.  Note  that  COMPLEX  type  is  not  present  in  ANSI-C 
standard.  It  is  defined  as  a  struct: 

typedef  struct  { 
double  real; 
double  imag; 

}  COMPLEX; 

and  thus  all  mathematical  operations  on  COMPLEX  type  variables  must  be  implemented  by 
the  user.  Because  the  model  describes  S-parameters  of  the  capacitor,  the  scattering  matrix  must 
be  converted  to  the  admittance  matrix  prior  return  statement.  The  conversion  is  done  by  the 
function  provided  by  ADS: 

extern  boolean  s_y_convert (COMPLEX  *inPar,  COMPLEX  *outPar,  int  direction,  double  rNo 

taking  as  the  parameters,  besides  input  and  output  matrices,  direction  of  conversion  (i.e.  S-to-Y 
or  Y-to-S),  reference  impedance  and  the  number  of  ports  of  the  model.  Complete  listing  of  the 
compute _y  function  of  MCM-D  capacitor  model  is  presented  below. 

static  boolean  compute_y( 

UserlnstDef  *userlnst  , 
double  omega, 

COMPLEX  *  yPar ) 

{ 

boolean  status  =  TRUE; 

COMPLEX  S matrix  [4]; 

double  f;  //  frequency 

double  x[6];  //  vector  of  parameters 

f  =  omega/2/PI;  //  conversion  to  Hz 

//  populating  ’  getmodelSXX  ’  input  vector 
x  [0]  =  f  / 1  e9  ; 

/* 

g_P  is  a  macro  : 

#define  G_P  userlnst  — >pData  [0]  .  value  .  dVal 
which  is  gap  width  parameter  in  milimeters 

*/ 
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x[l]  =  le3*g_P  ; 
x[2]  =  le3*w3_P; 
x[3]  =  le3*L_P; 
x[4]  =  le3*w2_P; 
x[5]  =  le3*wl_P; 

/* 

Function  : 

COMPLEX  getmodelSll(  int  n,  double*  x  ) 
return  sll  parameter  of  the  model,  where 
x  is  a  vector  of  parameters  in  meters  and  Hz 

*  / 

Smatrix[0]  =  Smatrix[3]  =  getmodelSll(  1,  &x[0]  ); 

Smatrix[l]  =  Smatrix[2]  =  getmodelS21(  1,  &x[0]  ); 

s_y_convert(  Smatrix  ,  yPar  ,  1,  ZO ,  2); 

return  status  ; 

} 

After  compiling  written  code  in  User-Compiled  Model  window  the  model  is  accessible  in 
Component  Library  in  schematic  window  (Fig. 2. 7). 


Figure  2.7:  View  of  Component  Library  window  in  schematic  showing  user-compiled  models. 


The  model  can  be  placed  into  schematic  window  in  the  same  way  as  built-in  models  and  used 
in  linear  simulation,  as  shown  on  Fig. 2. 8. The  results  of  simulation  of  the  circuit  in  schematic  as 
well  as  in  Momentum  (Fig. 2. 10)  are  shown  on  Fig. 2. 9. 
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Figure  2.8:  View  of  schematic  window  with  implemented  MCM-D  capacitor  model  and  ADS 
intergidital  capacitor  model. 


a)  b) 


Figure  2.9:  Results  of  simulations  of  the  example  of  MCM-D  capacitor.  The  red  line  represents 
ADS  interdigital  capacitor  model,  blue  line  represents  rational  model,  black  triangles  are  the 
results  of  Momentum  simulations. 

2.2.2  AWR  Microwave  Office 

Custom  models  in  Microwave  Office  can  be  created  using  Model  Wizard  from  Software  De¬ 
velopment  Kit.  The  modeling  wizard  generates  a  C++  description  of  the  model  that  can  be 
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Figure  2.10:  Layout  of  MCM-D  capacitor  in  Momentum. 


compiled  as  a  dynamically  linked  model.  The  models  are  implemented  using  Microsofts  COM 
(Component  Object  Model)  technology  used  for  communication  between  different  software  in 
Windows  family  operating  systems.  Conformity  to  COM  standard  allows  easy  migration  to 
future  versions  of  Microwave  Wizard  without  recompiling  any  binaries.  One  can  add  the  model 
just  by  placing  compiled  file  into  ’’Models”  subfolder  of  the  simulator  folder.  Steps  needed  to 
create  a  model  are  as  follows: 

1 .  Code  generation 

2.  Model  implementation 

3.  Compilation 

Similarly  to  Agilent’s  simulator,  most  of  the  source  code  is  generated  by  the  simulator,  precisely 
by  the  Model  Wizard  (see  Fig.??).  The  implementation  of  the  code  is  analogous  to  ADS  models 
and  is  limited  to  writing  the  mathematical  description  of  the  model  between  following  lines  of 
comments: 

// [-USER  CODE  BEGIN-] 

// [-USER  CODE  END-] 

Compilation  can  be  done  using  C++  compiler  linking  user-modified  .cpp  file  with  the  Compiler 
Dependent  Library.  Multiple  *.cpp  files  can  be  batch  compiled  so  that  multiple  models  can  be 
implemented  within  a  single  dll. 

Implementation  details 

KCmplxMtxData  Smat(2,2); 


double 

x  [6] ; 

x  [0]  = 

freq  / 1  e9 

x[l]  = 

g/le-3; 

x  [  2  ]  = 

w3  / 1  e  —  3; 

x[3]  = 

L/le-3; 

x  [4]  = 

w2/l  e— 3; 
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x  [5  ]  =  wl / 1  e  — 3; 


Smat(l,l)  =  Smat(2,2)  =  getmodelS  1 1  (  1,  &x[0]  ); 
Smat(l,2)  =  Smat(2,l)  =  getmodelS2 1  (  1,  &x[0]  ); 


AWR_Smat2Ymat  (  Smat ,  2,  NULL,  ymat  ); 


Figure  2.11:  Model  Wizard  in  Microwave  Office. 
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Figure  2.12:  Parameter  editor  in  the  Model  Wizard. 


Figure  2.13:  Model  Wizard  in  Microwave  Office. 
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Figure  2.14:  Model  Wizard  in  Microwave  Office. 
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